clc;
clear;
clf;

%% 区域几何尺寸及网格划分参数
H = 5; % 区域总高
L = 10; % 区域总长
Nx = 4; % 水平方向的网格数量
Ny = 2; % 竖直方向的网格数量, 选择能被 2 整除的数

%% 总单元数和节点总数
E = 2*Nx*Ny; % 总单元数
N = (Nx+1)*(Ny+1); % 节点总数

%% 单元间距
Dx = L/Nx; % 水平方向网格间距
Dy = H/Ny; % 竖直方向网格间距

%% 节点分布拓扑
AAA = zeros(Ny+1, Nx+1);
AAA(1,:) = 1:Nx+1;
for i = 2:Ny+1
    AAA(i,:) = AAA(i-1,:) + (Nx+1)*ones(1,Nx+1);
end

%% 三角形单元 JXY 生成
for i = 1:1*Ny+1
    for j = 1:1*Nx+1
        JXY(AAA(i,j), 1) = Dx*(j-1);
        JXY(AAA(i,j), 2) = Dy*(i-1);
    end
end

%% 网格平面旋转
t = 0;
for i = 1:length(JXY(:,1))
    R = sqrt((JXY(i,1)+1)^2 + JXY(i,2)^2);
    t1 = atan(JXY(i,2)/(JXY(i,1)+1));
    JXY(i,1) = R*cos(t/180*pi + t1);
    JXY(i,2) = R*sin(t/180*pi + t1);
end

%% 三角形单元 JM 生成
k = 0;
for i = 1:Ny/2
    for j = 1:Nx
        k = k + 1;
        JM(k,:) = [AAA(i,j), AAA(i,j+1), AAA(i+1,j)];
        k = k + 1;
        JM(k,:) = [AAA(i+1,j+1), AAA(i+1,j), AAA(i,j+1)];
    end
end
for i = Ny/2+1:Ny
    for j = 1:Nx
        k = k + 1;
        JM(k,:) = [AAA(i+1,j+1), AAA(i,j), AAA(i,j+1)];
        k = k + 1;
        JM(k,:) = [AAA(i,j), AAA(i+1,j+1), AAA(i+1,j)];
    end
end

%% BP 数据生成
BP1 = AAA(1,:);
BP2 = AAA(:,Nx+1);
BP3 = AAA(Ny+1,:);
BP4 = AAA(:,1);

%% BE 数据生成
tx1 = pi/2 - t/180*pi; % 1 号边界外法线方向与 x 轴夹角
ty1 = pi - t/180*pi; % 1 号边界外法线方向与 y 轴夹角
tx2 = t/180*pi; % 2 号边界外法线方向与 x 轴夹角
ty2 = pi/2 - tx2; % 2 号边界外法线方向与 y 轴夹角
tx3 = pi - pi/2 + t/180*pi; % 3 号边界外法线方向与 x 轴夹角
ty3 = pi - t/180*pi + pi; % 3 号边界外法线方向与 y 轴夹角
tx4 = (180 + t)/180*pi; % 4 号边界外法线方向与 x 轴夹角
ty4 = pi/2 + t/180*pi; % 4 号边界外法线方向与 y 轴夹角
BE1 = [(1:2:2*Nx)', ones(size((1:Nx)')), ones(Nx,1)*cos(tx1), ones(Nx,1)*cos(ty1)];
BBB = [(2*Nx:2*Nx:Ny*Nx)'; (Ny*Nx+2*Nx-1:2*Nx:2*Ny*Nx)'];
BE2 = [BBB, 3*ones(size((1:Ny)')), ones(Ny,1)*cos(tx2), ones(Ny,1)*cos(ty2)];
CCC = (2*Nx*(Ny-1)+2:2:2*Nx*Ny)';
BE3 = [CCC, 2*ones(size((1:Nx)')), ones(Nx,1)*cos(tx3), ones(Nx,1)*cos(ty3)];
DDD = [(1:2*Nx:Nx*(Ny-1)+1)'; (Nx*Ny+2:2*Nx:2*Nx*Ny)'];
BE4 = [DDD, 3*ones(size((1:Ny)')), ones(Ny,1)*cos(tx4), ones(Ny,1)*cos(ty4)];

%% 调用三角形网格绘制程序
triangle_grid(JM, JXY);

%% 清除多余变量, 存储网格数据
clear Dx Dy
clear H L Nx Ny i j k
clear t t1 R
clear AAA BBB CCC DDD
clear tx1 tx2 tx3 tx4
clear ty1 ty2 ty3 ty4
save msh

